Fan-beam and cone-beam image reconstruction using filtered backprojection of differentiated projection data

ABSTRACT

A tomographic image reconstruction method produces either 2D or 3D images from fan beam or cone beam projection data by filtering the backprojection image of differentiated projection data. The reconstruction is mathematically exact if sufficient projection data is acquired. A cone beam embodiment and both a symmetric and asymmetric fan beam embodiment are described.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is based on U.S. Provisional Patent Application Ser.No. 60/630,930, filed on Nov. 24, 2004 and entitled “FAN-BEAM ANDCONE-BEAM IMAGE RECONSTRUCTION USING FILTERED BACKPROJECTION OFDIFFERENTIATED PROJECTION DATA,” and U.S. Provisional Patent ApplicationSer. No. 60/630,932, filed on Nov. 24, 2004 and entitled “FAN-BEAM IMAGERECONSTRUCTION METHOD FOR TRUNCATED PROJECTION DATA.”

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

This invention was made with government support under Grant No. EB001683awarded by the National Institute of Health. The United StatesGovernment has certain rights in this invention.

BACKGROUND OF THE INVENTION

The present invention relates to computed tomography (CT) imagingapparatus; and more particularly, to a method for reconstructing imagesfrom divergent beams of acquired image data.

In a current computed tomography system, an x-ray source projects afan-shaped beam which is collimated to lie within an X-Y plane of aCartesian coordinate system, termed the “imaging plane.” The x-ray beampasses through the object being imaged, such as a medical patient, andimpinges upon an array of radiation detectors. The intensity of thetransmitted radiation is dependent upon the attenuation of the x-raybeam by the object and each detector produces a separate electricalsignal that is a measurement of the beam attenuation. The attenuationmeasurements from all the detectors are acquired separately to producethe transmission profile.

The source and detector array in a conventional CT system are rotated ona gantry within the imaging plane and around the object so that theangle at which the x-ray beam intersects the object constantly changes.A group of x-ray attenuation measurements from the detector array at agiven angle is referred to as a “view” and a “scan” of the objectcomprises a set of views made at different angular orientations duringone revolution of the x-ray source and detector. In a 2D scan, data isprocessed to construct an image that corresponds to a two dimensionalslice taken through the object. The prevailing method for reconstructingan image from 2D data is referred to in the art as the filteredbackprojection technique. This process converts the attenuationmeasurements from a scan into integers called “CT numbers” or“Hounsfield units”, which are used to control the brightness of acorresponding pixel on a display.

The term “generation” is used in CT to describe successivelycommercially available types of CT systems utilizing different modes ofscanning motion and x-ray detection. More specifically, each generationis characterized by a particular geometry of scanning motion, scanningtime, shape of the x-ray beam, and detector system.

As shown in FIG. 1, the first generation utilized a single pencil x-raybeam and a single scintillation crystal-photomultiplier tube detectorfor each tomographic slice. After a single linear motion or traversal ofthe x-ray tube and detector, during which time 160 separate x-rayattenuation or detector readings are typically taken, the x-ray tube anddetector are rotated through 1° and another linear scan is performed toacquire another view. This is repeated typically to acquire 180 views.

A second generation of devices developed to shorten the scanning timesby gathering data more quickly is shown in FIG. 2. In these units amodified fan beam in which anywhere from three to 52 individualcollimated x-ray beams and an equal number of detectors are used.Individual beams resemble the single beam of a first generation scanner.However, a collection of from three to 52 of these beams contiguous toone another allows multiple adjacent cores of tissue to be examinedsimultaneously. The configuration of these contiguous cores of tissueresembles a fan, with the thickness of the fan material determined bythe collimation of the beam and in turn determining the slice thickness.Because of the angular difference of each beam relative to the others,several different angular views through the body slice are beingexamined simultaneously. Superimposed on this is a linear translation orscan of the x-ray tube and detectors through the body slice. Thus, atthe end of a single translational scan, during which time 160 readingsmay be made by each detector, the total number of readings obtained isequal to the number of detectors times 160. The increment of angularrotation between views can be significantly larger than with a firstgeneration unit, up to as much as 36°. Thus, the number of distinctrotations of the scanning apparatus can be significantly reduced, with acoincidental reduction in scanning time. By gathering more data pertranslation, fewer translations are needed.

To obtain even faster scanning times it is necessary to eliminate thecomplex translational-rotational motion of the first two generations. Asshown in FIG. 3, third generation scanners therefore use a much widerfan beam. In fact, the angle of the beam may be wide enough to encompassmost or all of an entire patient section without the need for a lineartranslation of the x-ray tube and detectors. As in the first twogenerations, the detectors, now in the form of a large array, arerigidly aligned relative to the x-ray beam, and there are notranslational motions at all. The tube and detector array aresynchronously rotated about the patient through an angle of 180-360°.Thus, there is only one type of motion, allowing a much faster scanningtime to be achieved. After one rotation, a single tomographic section isobtained.

Fourth generation scanners feature a wide fan beam similar to the thirdgeneration CT system as shown in FIG. 4. As before, the x-ray tuberotates through 360° without having to make any translational motion.However, unlike in the other scanners, the detectors are not alignedrigidly relative to the x-ray beam. In this system only the x-ray tuberotates. A large ring of detectors are fixed in an outer circle in thescanning plane. The necessity of rotating only the tube, but not thedetectors, allows faster scan time.

Most of the commercially available CT systems employ imagereconstruction methods based on the concepts of Radon space and theRadon transform. For the pencil beam case, the data is automaticallyacquired in Radon space. Therefore a Fourier transform can directlysolve the image reconstruction problem by employing the well-knownFourier-slice theorem. Such an image reconstruction procedure is calledfiltered backprojection (FBP). The success of FBP reconstruction is dueto the translational and rotational symmetry of the acquired projectiondata. In other words, in a parallel beam data acquisition, theprojection data are invariant under a translation and/or a rotationabout the object to be imaged. For the fan beam case, one can solve theimage reconstruction problem in a similar fashion, however, to do thisan additional “rebinning” step is required to transform the fan beamdata into parallel beam data. The overwhelming acceptance of theconcepts of Radon space and the Radon transform in the two dimensionalcase gives this approach to CT image reconstruction a paramount positionin tomographic image reconstruction.

The Radon space and Radon transformation reconstruction methodology ismore problematic when applied to three-dimensional image reconstruction.Three-dimensional CT, or volume CT, employs an x-ray source thatprojects a cone beam on a two-dimensional array of detector elements asshown in FIG. 5. Each view is thus a 2D array of x-ray attenuationmeasurements and a complete scan produced by acquiring multiple views asthe x-ray source and detector array are revolved around the subjectresults in a 3D array of attenuation measurements. The reason for thisdifficulty is that the simple relation between the Radon transform andthe x-ray projection transform for the 2D case in not valid in the 3Dcone beam case. In the three-dimensional case, the Radon transform isdefined as an integral over a plane, not an integral along a straightline. Consequently, we have difficulty generalizing the success of theRadon transform as applied to the 2D fan beam reconstruction to the 3Dcone beam reconstruction. In other words, we have not managed to derivea shift-invariant FBP method by directly rebinning the measured conebeam data into Radon space. Numerous solutions to this problem have beenproposed as exemplified in U.S. Pat. Nos. 5,270,926; 6,104,775;5,257,183; 5,625,660; 6,097,784; 6,219,441; and 5,400,255.

It is well known that the projection-slice theorem (PST) plays animportant role in the image reconstruction from two- andthree-dimensional parallel-beam projections. The power of the PST liesin the fact that Fourier transform of a single view of parallel-beamprojections is mapped into a single line (two-dimensional case) or asingle slice (three-dimensional case) in the Fourier space via the PST.In other words, a complete Fourier space of the image object can bebuilt up from the Fourier transforms of the sequentially measuredparallel-beam projection data. Once all the Fourier information of theimage object is known, an inverse Fourier transform can be performed toreconstruct the image. Along the direction of the parallel-beamprojections, there is a shift-invariance of the image object in a singleview of the parallel-beam projections. This is the fundamental reasonfor the one-to-one correspondence between the Fourier transform ofparallel-beam projections and a straight line or a slice in the Fourierspace. The name of the projection-slice theorem follows from thisone-to-one correspondence.

In practice, divergent fan-beam and cone-beam scanning modes have thepotential to allow fast data acquisition. But image reconstruction fromdivergent-beam projections poses a challenge. In particular, the PST isnot directly applicable to the divergent-beam projections since theshift-invariance in a single view of projections is lost in thedivergent-beam cases. One way to bypass this problem is to explicitlyrebin the measured divergent-beam projections into parallel beamprojections. This is the basic method currently used in solving theproblem of fan-beam image reconstruction. After the rebinning process,one can take the advantages of the fast Fourier transforms (FFT) toefficiently reconstruct images. There are some issues on the potentialloss of image spatial resolution due to the data rebinning. But thereare also some advantages in generating uniform distribution of imagenoise due to the non-local characteristic of the Fourier transform.Alternatively, a fan-beam projection can also be relabeled in terms ofRadon variables so that the two-dimensional inverse Radon transform canbe used to reconstruct images. In this way, a convolution-based fan-beamimage reconstruction algorithm can be readily developed. The advantageof this type of reconstruction algorithm is the explicit filteredbackprojection (FBP) structure. The disadvantage of theconvolution-based method is that the weight in the backprojection stepdepends on the individual image pixels and thus noise distribution maynot be uniform. This may pose problems in the clinical interpretation oftomographic images. In practice, different CT manufactures may utilizedifferent strategies in balancing these advantages and disadvantages.

In the cone-beam case, it is much more complicated to rebin cone-beamprojections into parallel-beam projections. The huge cone-beam data setalso poses a big challenge to the potential data storage during therebinning process. The main stream of the developments in cone-beamreconstruction has been focused on the development of approximate orexact reconstruction methods. For circular-based source trajectories,methods disclosed by L. A. Feldkamp, L. C. Davis, and J. W. Kress,“Practical Cone Beam Algorithm,” J. Opt. Soc. Am. A 1, 612-619(1984); G.Wang, T. H. Lin, P. Cheng, and D. M. Shinozaki, “A general cone-beamreconstruction algorithm,” IEEE Trans. Med. Imaging 12, 486-496 (1993);generate acceptable image quality up to moderate cone angles (up to 10°or so). Exact reconstruction algorithms have also been proposed andfurther developed for both helical source trajectory and more generalsource trajectories. Most recently, a mathematically exact andshift-invariant FBP reconstruction formula was proposed for thehelical/spiral source trajectory A. Katsevich, “Theoretically exactfiltered backprojection-type inversion algorithm for spiral CT,” SIAM(Soc. Ind. Appl. Math.) J. Appl. Math. 62, 2012-2026 (2002). Startingwith either the original Tuy's framework or Grangeat's framework, uponan appropriate choice of weighting function over the redundant data,shift-invariant FBP reconstruction formula has been derived for ageneral source trajectory. Similar to the fan-beam FBP reconstructionalgorithm the characteristic of the convolution-based cone-beamreconstruction algorithm is the voxel-dependent weighting factor in thebackprojection step. This will cause non-uniform distribution of theimage noise. Moreover, due to the local nature of the newly developedconvolution-based cone-beam image reconstruction algorithms, differentimage voxels are reconstructed by using cone-beam projection dataacquired at different pieces of the source trajectory. Namely, differentimage voxels are reconstructed by using the data acquired underdifferent physical conditions. This will potentially lead to some datainconsistency in dynamic imaging. Finally, the current convolution-basedimage reconstruction algorithms are only valid for some discrete pitchvalues in the case of helical/spiral source trajectory. This featurelimits their application in a helical/spiral cone-beam CT scanner.

SUMMARY OF THE INVENTION

The present invention is an image reconstruction method which can begenerically defined as an integral of a product of two functions,K({right arrow over (x)},{right arrow over (x)}′) and Q({right arrowover (x)},{right arrow over (x)}′). The function Q({right arrow over(x)},{right arrow over (x)}′) is a weighted backprojection operation onthe differentiated projection data acquired during a scan. The functionK({right arrow over (x)},{right arrow over (x)}′) is the Fouriertransform of one of the weighting factors, ω₁({right arrow over(x)},{circumflex over (k)}), of the factorized weighting functionω({right arrow over (x)},{circumflex over (k)};t). An efficientreconstruction method results when the second weighting factor ω₂({rightarrow over (x)},t) of the weighting function ({right arrow over(x)},{circumflex over (k)}; t) is chosen to be a piecewise constant, orin other words, a function consisting of several constant pieces. Inthis case, the filtering kernel K({right arrow over (x)},{right arrowover (x)}′) is reduced to a Hilbert kernel or a linear combination ofHilbert kernels along different filtering lines.

More specifically, the invention is a method for producing an image froma plurality of projection views of a subject acquired at different viewangles which includes: differentiating the acquired projection views;producing image values along filter lines in image space bybackprojecting differentiated projection views; and filtering imagevalues along each filter line. The invention may be employed with manydifferent imagining modalities that use either symmetric or asymmetricdiverging beams and either 2D fan beams or 3D cone beams.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1-4 are schematic representations of different CT systemgeometries;

FIG. 5 is a pictorial representation of a 3D, or volume, CT system;

FIG. 6 is a pictorial view of a CT system which employs the presentinvention;

FIG. 7 is a block diagram of the CT system of FIG. 6;

FIGS. 8 a, 8 b, 9, 10 and 11 are pictorial representations of the sourcetrajectory and the field of view of the reconstructed image which areused to illustrate the underlying theory of the present invention;

FIG. 12 is a flow chart of the steps performed by the imaging system ofFIGS. 6 and 7 to practice a preferred embodiment of the presentinvention; and

FIG. 13 is a pictorial representation of an asymmetric detector used inan alternative embodiment of the invention.

GENERAL DESCRIPTION OF THE INVENTION

A compactly supported function f({right arrow over (x)}) is the imagefunction to be reconstructed. A general source trajectory isparameterized by a single parameter t, t∈Λ=[t_(i), t_(f)]. Here t_(i)and t_(f) are assumed to be the view angles corresponding to thestarting and ending points of an open source trajectory. A point on thesource trajectory is measured in a laboratory system and denoted as{right arrow over (y)}(t). The divergent fan-beam and cone-beamprojections from a source point {right arrow over (y)}(t) are definedas:

$\begin{matrix}{{g\left\lbrack {\overset{->}{r},{\overset{->}{y}(t)}} \right\rbrack} = {\int_{0}^{+ \infty}{{\mathbb{d}{{sf}\left\lbrack {{\overset{->}{y}(t)} + {s\overset{->}{r}}} \right\rbrack}}.}}} & (2.1)\end{matrix}$

An intermediate function G_(D)[{right arrow over (k)},{right arrow over(y)}(t)] is defined as the Fourier transform of the above divergent beamprojections g[{right arrow over (r)},{right arrow over (y)}(t)] withrespect to the variable {right arrow over (r)}, viz.,

$\begin{matrix}{{G_{D}\left\lbrack {\overset{->}{k},{\overset{->}{y}(t)}} \right\rbrack} = {\int_{R^{D}}{{\mathbb{d}\overset{->}{r}}{g\left\lbrack {\overset{->}{r},{\overset{->}{y}(t)}} \right\rbrack}{{\mathbb{e}}^{{\mathbb{i}}\; 2\;\pi\;{\overset{->}{k} \cdot \overset{->}{r}}}.}}}} & (2.2)\end{matrix}$The letter D labels the dimensions of a Euclidean space. Thus, D=2corresponds to the fan-beam case, and D=3 corresponds to the cone-beamcase.

The convention of decomposing a vector into its magnitude and a unitdirection vector will be used, i.e. {right arrow over (r)}=r{circumflexover (r)} and {right arrow over (k)}=k{circumflex over (k)}. Using theabove two definitions, the following properties of the projections andintermediate functions can be readily demonstrated:

$\begin{matrix}{{{g\left\lbrack {\overset{\rightarrow}{r},{\overset{\rightarrow}{y}(t)}} \right\rbrack} = {\frac{1}{r}{\overset{\_}{g}\left\lbrack {\hat{r},{\overset{\rightarrow}{y}(t)}} \right\rbrack}}},} & (2.3) \\{{{G_{D}\left\lbrack {\overset{\rightarrow}{k},{\overset{\rightarrow}{y}(t)}} \right\rbrack} = {\frac{1}{\kappa^{D - 1}}{{\overset{\_}{G}}_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(t)}} \right\rbrack}}},} & (2.4) \\{{{{Im}\;{G_{D}\left\lbrack {{- \overset{\rightarrow}{k}},{\overset{\rightarrow}{y}(t)}} \right\rbrack}} = {{- {Im}}\;{G_{D}\left\lbrack {\overset{\rightarrow}{k},{\overset{\rightarrow}{y}(t)}} \right\rbrack}}},} & (2.5)\end{matrix}$where Im G_(D)[{right arrow over (k)},{right arrow over (y)}(t)] is theimaginary part of the intermediate function G_(D)[−{right arrow over(k)},{right arrow over (y)}(t)]. In addition, g[{circumflex over(r)},{right arrow over (y)}(t)] and G _(D)[{circumflex over (k)},{rightarrow over (y)}(t)] are the angular components of the projection g[{right arrow over (r)},{right arrow over (y)}(t)] and the intermediatefunction G_(D)[{right arrow over (k)},{right arrow over (y)}(t)]respectively.

With the above definitions and mathematical convention in hand, theimage function f({right arrow over (x)}) may be reconstructed using themodified Tuy-like formula:

$\begin{matrix}{{f\left( \overset{\rightarrow}{x} \right)} = {{\int_{S_{\kappa}^{D - 1}}^{\mspace{11mu}}\ {{\mathbb{d}\hat{k}}{\sum\limits_{j}{\frac{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t_{j}}} \right)}{2\;\pi\;{\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}\left( t_{j} \right)}}}\frac{\partial\;}{\partial_{q}}{Im}\;{{\overset{\_}{G}}_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}}}}❘_{q = t_{j}}.}} & (2.6)\end{matrix}$Here ({right arrow over (x)},{circumflex over (k)}; t_(j)) is aweighting function used in order to take into account the multiplesolutions of the following equation:{circumflex over (k)}·[{right arrow over (x)}−{right arrow over (y)}(t_(j))]=0  (2.7)In the cone-beam case, the unit vector {circumflex over (k)} is thenormal vector of a plane that contains the line connecting a given imagepoint ({right arrow over (x)}), and a source point ({right arrow over(y)}(t_(j))). While, in the fan-beam case the unit vector {circumflexover (k)} is the normal vector to the line connecting a given imagepoint ({right arrow over (x)}), and a source point ({right arrow over(y)}(t_(j))). This equation dictates the well-known Tuy data sufficiencycondition: Any plane (cone-beam case) or straight line (fan-beam case)that passes through the region of interest (ROI) should intersect thesource trajectory at least once. In Eq. 2.6 and Eq. 2.7, {right arrowover (y)}(t_(j)) labels the j-th intersection point between the planecontaining the image point ({right arrow over (x)}) and the sourcetrajectory {right arrow over (y)}(t_(j)). The set which is the union ofall these solutions is labeled as T({right arrow over (x)},{circumflexover (k)}), viz.,t _(j) ∈T({right arrow over (x)},{right arrow over (k)})={t|{circumflexover (k)}·[{right arrow over (x)}−{right arrow over (y)}(t)]=0}.  (2.8)The weighting function ω({right arrow over (x)},{circumflex over (k)};t) satisfies the following normalization condition:

$\begin{matrix}{{\sum\limits_{t_{j} \in {({\overset{\rightarrow}{x},\hat{k}})}}{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t_{j}}} \right)}} = 1.} & (2.9)\end{matrix}$In addition, using property (2.5) and the inversion formula (2.6), theweighting function possesses the following symmetry:ω({right arrow over (x)},−{circumflex over (k)};t)=ω({right arrow over(x)},{circumflex over (k)};t)  (2.10)In order to proceed, we use the following well-known formula

$\begin{matrix}{{\sum\limits_{t_{j}}\frac{\delta\left( {t - t_{j}} \right)}{{\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}\left( t_{j} \right)}}}} = {\delta\left\{ {\hat{k} \cdot \left\lbrack {\overset{\rightarrow}{x} - {\overset{\rightarrow}{y}(t)}} \right\rbrack} \right\}}} & (2.11)\end{matrix}$to turn the summation over t_(j) into an integral over the continuousvariable t in Eq. (2.6). Consequently, we obtain:

$\begin{matrix}{{f\left( \overset{\rightarrow}{x} \right)} = {{\frac{1}{2\;\pi}{\int_{\Lambda}^{\;}\mspace{7mu}{{\mathbb{d}t}{\int_{S_{\kappa}^{D - 1}}^{\;}\mspace{7mu}{{\mathbb{d}\hat{k}}\;{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t}} \right)}{{sgn}\left\lbrack {\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}(t)}} \right\rbrack}\frac{\partial\;}{\partial q}{Im}\;{{\overset{\_}{G}}_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}}}}}❘_{q = t}{{\delta\left\lbrack {\hat{k} \cdot \left( {\overset{\rightarrow}{x} - {\overset{\rightarrow}{y}(t)}} \right)} \right\rbrack}.}}} & (2.12)\end{matrix}$Starting with the above equation, one can derive a shift-invariant FBPreconstruction formula for the fan-beam and cone-beam projectionsrespectively. However, in order to derive the desired FBPD formula, itis insightful to use the following integration representation of thedelta function:

$\begin{matrix}{{\delta\left\lbrack {\overset{\bigwedge}{k} \cdot \left( {\overset{->}{x} - {\overset{->}{y}(t)}} \right)} \right\rbrack} = {\int_{- \infty}^{+ \infty}{{\mathbb{d}k}\;{\mathbb{e}}^{{\mathbb{i}}\; 2\;\pi\; k{\hat{k} \cdot {({\overset{->}{x} - {\overset{->}{y}{(t)}}})}}}}}} & (2.13)\end{matrix}$to rewrite Eq. (2.12) as follows

$\begin{matrix}{{f\left( \overset{\rightarrow}{x} \right)} = {{\frac{1}{2\;\pi}{\int_{\Lambda}^{\;}\mspace{7mu}{{\mathbb{d}t}{\int_{0}^{+ \infty}\mspace{7mu}{{\mathbb{d}k}{\int_{S^{D - 1}}^{\;}\mspace{7mu}{{\mathbb{d}\hat{k}}{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t}} \right)}{{sgn}\left\lbrack {\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}(t)}} \right\rbrack}\frac{\partial\;}{\partial q}{Im}\;{{\overset{\_}{G}}_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}}}}}}}❘_{q = t}{{\mathbb{e}}^{{\mathbb{i}2}\;\pi\; k{\hat{k} \cdot {\lbrack{\overset{\rightarrow}{x} - {\overset{\rightarrow}{y}{(t)}}}\rbrack}}} + {\frac{1}{2\;\pi}{\int_{\Lambda}^{\;}\mspace{7mu}{{\mathbb{d}t}{\int_{- \infty}^{0}\mspace{7mu}{{\mathbb{d}k}{\int_{S^{D - 1}}^{\;}\ {{\mathbb{d}\hat{k}}{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t}} \right)}{sgn}{\quad{{\left\lbrack {\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}(t)}} \right\rbrack\frac{\partial\;}{\partial q}{Im}\;{{\overset{\_}{G}}_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}❘_{q = t}{{\mathbb{e}}^{{\mathbb{i}2}\;\pi\; k{\hat{k} \cdot {\lbrack{\overset{\rightarrow}{x} - {\overset{\rightarrow}{y}{(t)}}}\rbrack}}}.}}}}}}}}}}}}} & (2.14)\end{matrix}$By changing the variables k→−k and {circumflex over (k)}→−{circumflexover (k)} in the second term and using Eqs. (2.10) and (2.5), we obtain:

$\begin{matrix}{{f\left( \overset{\rightarrow}{x} \right)} = {\frac{1}{\;\pi}{\int_{\Lambda}^{\;}\mspace{7mu}{{\mathbb{d}t}{\int_{0}^{+ \infty}\mspace{7mu}{{\mathbb{d}k}{\int_{S^{D - 1}}^{\;}\mspace{7mu}{{\mathbb{d}\hat{k}}{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t}} \right)}{{sgn}\left\lbrack {\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}(t)}} \right\rbrack}\frac{\partial\;}{\partial q}{Im}\;{\overset{\_}{G}}_{D}{\quad{\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack ❘_{q = t}{\mathbb{e}}^{{\mathbb{i}2}\;\pi\; k{\hat{k} \cdot {\lbrack{\overset{\rightarrow}{x} - {\overset{\rightarrow}{y}{(t)}}}\rbrack}}}}}}}}}}}}} & (2.15)\end{matrix}$In this equation, the dummy variable k takes only non-negative values.Thus, the unit vector {circumflex over (k)} and the non-negative value kcan be combined as a complete vector, i.e., {right arrow over(k)}=k{circumflex over (k)}. Multiplying and dividing the angularcomponent Im G_(D)[{circumflex over (k)},{right arrow over (y)}(t)] by afactor k^(D-1) and using property (2.4), we obtain:

$\begin{matrix}{{{f\left( \overset{\rightarrow}{x} \right)} = {{\frac{1}{\;\pi}{\int_{\Lambda}^{\;}\mspace{7mu}{{\mathbb{d}t}{\int_{R_{k}^{D}}^{\;}\mspace{7mu}{{\mathbb{d}\overset{\rightarrow}{k}}\;{\omega\left( {\overset{\rightarrow}{x},{\hat{k};t}} \right)}{{sgn}\left\lbrack {\hat{k} \cdot {{\overset{\rightarrow}{y}}^{\prime}(t)}} \right\rbrack}\frac{\partial\;}{\partial q}{Im}\;{G_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}}}}}❘_{q = t}{\mathbb{e}}^{{\mathbb{i}2}\;\pi\; k{\hat{k} \cdot {\lbrack{\overset{\rightarrow}{x} - {\overset{\rightarrow}{y}{(t)}}}\rbrack}}}}},} & (2.16)\end{matrix}$where we have used the volume element d{right arrow over(k)}=dkk^(D-1)d{circumflex over (k)} in a D-dimensional Euclidean space.

In the following, we use the definition of the intermediate functionG_(D)[{right arrow over (k)},{right arrow over (y)}(t)] to calculate thederivative with respect to the source parameter:

$\begin{matrix}{{{\frac{\partial\;}{\partial q}{Im}\;{G_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}❘_{q = t}} = {{{Im}{\int_{R^{D}}^{\;}\mspace{7mu}{{\mathbb{d}\overset{\rightarrow}{r}}\frac{\partial\;}{\partial q}{g\left\lbrack {\overset{\rightarrow}{r},{\overset{\rightarrow}{y}(q)}} \right\rbrack}_{q = t}{\mathbb{e}}^{{- {\mathbb{i}2}}\;\pi\;{\overset{\rightarrow}{k} \cdot \overset{\rightarrow}{r}}}}}} = {{\int_{R^{D}}^{\;}\mspace{7mu}{{\mathbb{d}\overset{\rightarrow}{r}}\frac{\partial\;}{\partial q}{g\left\lbrack {\overset{\rightarrow}{r},{\overset{\rightarrow}{y}(q)}} \right\rbrack}_{q = t}{{Im}\left\lbrack {\mathbb{e}}^{{\mathbb{i}2}\;\pi\;{\overset{\rightarrow}{k} \cdot \overset{\rightarrow}{r}}} \right\rbrack}}} = {{\int_{R^{D}}^{\;}\mspace{7mu}{{\mathbb{d}\overset{\rightarrow}{r}}\frac{\partial\;}{\partial q}{g\left\lbrack {\overset{\rightarrow}{r},{\overset{\rightarrow}{y}(q)}} \right\rbrack}_{q = t}\frac{1}{2{\mathbb{i}}}\left( {{\mathbb{e}}^{{- {\mathbb{i}2}}\;\pi\;{\overset{\rightarrow}{k} \cdot \overset{\rightarrow}{r}}} - {\mathbb{e}}^{{- {\mathbb{i}2}}\;\pi\;{\overset{\rightarrow}{k} \cdot \overset{\rightarrow}{r}}}} \right)}} = {{\frac{1}{2{\mathbb{i}}}{\int_{R^{D}}^{\;}\mspace{7mu}{{\mathbb{d}\overset{\rightarrow}{r}}\frac{\partial\;}{\partial q}\left\{ {{g\left\lbrack {\overset{\rightarrow}{r},{\overset{\rightarrow}{y}(q)}} \right\rbrack} - {g\left\lbrack {{- \overset{\rightarrow}{r}},{\overset{\rightarrow}{y}(q)}} \right\rbrack}} \right\}}}}❘_{q = t}{{\mathbb{e}}^{{- {\mathbb{i}2}}\;\pi\;{\overset{\rightarrow}{k} \cdot \overset{\rightarrow}{r}}}.}}}}}} & (2.17)\end{matrix}$Note that the variable {right arrow over (r)} is defined in a localcoordinate system, i.e.,{right arrow over (r)}={right arrow over (x)}′−{right arrow over(y)}(t),  (2.18)where {right arrow over (x)}′ is a point in the backprojection imagespace. Therefore, for the fan-beam and cone-beam projections at a givensource position {right arrow over (y)}(t), the substitution of Eq.(2.18) into Eq. (2.17), and the application of property (2.3) yields:

$\begin{matrix}{{{{\frac{\partial\;}{\partial q}{Im}\;{{\overset{\_}{G}}_{D}\left\lbrack {\hat{k},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}❘_{q = t}} = {{\frac{1}{2{\mathbb{i}}}{\int_{R^{D}}^{\;}\mspace{7mu}{{\mathbb{d}{\overset{\rightarrow}{x}}^{\prime}}\frac{1}{{{\overset{\rightarrow}{x}}^{\prime} - {\overset{\rightarrow}{y}(t)}}}\frac{\partial\;}{\partial q}{{\overset{\_}{g}}^{A}\left\lbrack {{\hat{\beta}}_{\overset{\rightarrow}{x}},{\overset{\rightarrow}{y}(q)}} \right\rbrack}}}}❘_{q = t}{\mathbb{e}}^{{- {\mathbb{i}}}\; 2\pi{\overset{\rightarrow}{k} \cdot {\lbrack{{\overset{\rightarrow}{x}}^{\prime} - {\overset{\rightarrow}{y}{(t)}}}\rbrack}}}}},} & (2.19)\end{matrix}$where the unit vector {circumflex over (β)}_({right arrow over (x)}), isdefined as:

$\begin{matrix}{{\hat{\beta}}_{{\overset{->}{x}}^{\prime}} = {\frac{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}{{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}}.}} & (2.20)\end{matrix}$The modified fan-beam and cone-beam data g ^(A)[{circumflex over(β)}_({right arrow over (x)}′),{right arrow over (y)}(t)] is defined as:g ^(A)[{circumflex over (β)}_({right arrow over (x)}′) ,{right arrowover (y)}(t)]= g[{circumflex over (β)} _({right arrow over (x)}′),{right arrow over (y)}(t)]− g[−{circumflex over (β)}_({right arrow over (x)}′) ,{right arrow over (y)}(t)].  (2.21)The meaning of g ^(A)[{circumflex over(β)}_({right arrow over (x)}′),{right arrow over (y)}(t)] will beelaborated upon later. After Eq. (2.16) and switching the order ofintegrations, we obtain:

$\begin{matrix}{{f\mspace{11mu}\left( \overset{->}{x} \right)} = {{\frac{1}{2{\pi\mathbb{i}}}{\int_{R^{D}}\ {{\mathbb{d}{\overset{->}{x}}^{\prime}}{\int_{R_{k}^{D}}\ {{\mathbb{d}\overset{->}{k}}{\mathbb{e}}^{{\mathbb{i}2\pi}{\overset{->}{k} \cdot {({\overset{->}{x} - {\overset{->}{x}}^{\prime}})}}}{\int_{\Lambda}{{\mathbb{d}t}\frac{\omega\mspace{11mu}\left( {\overset{->}{x},{\hat{k};t}} \right)\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot {{\overset{->}{y}}^{\prime}(t)}} \right\rbrack}}{{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}}\frac{\partial}{\partial q}{{\overset{\_}{g}}^{A}\left\lbrack {{\hat{\beta}}_{{\overset{->}{x}}^{\prime}},{\overset{->}{y}\mspace{11mu}(q)}} \right\rbrack}}}}}}}}❘_{q = t}.}} & (2.22)\end{matrix}$In the above equation, the integration over the variables t and {rightarrow over (k)} are coupled to one another through the factor ω({rightarrow over (x)},{circumflex over (k)};t)sgn[{circumflex over (k)}·{rightarrow over (y)}′(t)]. A key observation is to impose the followingconstraint, i.e., a factorization property, on the choice of theweighting functions:ω({right arrow over (x)},{circumflex over (k)};t)=ω₁({right arrow over(x)},{circumflex over (k)})ω₂({right arrow over (x)},t)sgn[{circumflexover (k)}·{square root over (y)}′(t)].  (2.23)In other words, we require that the weighting function is factorable sothat the integrations over the variables t and {right arrow over (k)}can be decoupled. This factorized weighting function is an importantaspect of the present invention. Under this factorization condition, thereconstruction formula (2.22) can be further simplified as:

$\begin{matrix}{{{f\mspace{11mu}\left( \overset{->}{x} \right)} = {\int_{R^{D}}\ {{\mathbb{d}{\overset{->}{x}}^{\prime}}K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)\mspace{11mu} Q\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)}}}{where}} & (2.24) \\{{{K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {\frac{1}{2{\pi\mathbb{i}}}{\int_{R_{k}^{D}}\ {{\mathbb{d}\overset{->}{k}}{\omega_{1}\left( {\overset{->}{x},\hat{k}} \right)}\mspace{11mu}{\mathbb{e}}^{{\mathbb{i}2\pi}{\overset{->}{k} \cdot {({\overset{->}{x} - {\overset{->}{x}}^{\prime}})}}}}}}},} & (2.25) \\{{{Q\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = \;{{\int_{\Lambda}{{\mathbb{d}t}\frac{\omega_{2}\;\left( {\overset{->}{x},t} \right)}{{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}}\frac{\partial}{\partial q}{{\overset{\_}{g}}^{A}\left\lbrack {{\hat{\beta}}_{{\overset{->}{x}}^{\prime}},{\overset{->}{y}\mspace{11mu}(q)}} \right\rbrack}}}❘_{q = t}.}}\;} & (2.26)\end{matrix}$

In Eq. (2.24), the image reconstruction consists of two major steps. Inthe first step, the fan-beam or cone-beam projections are weighted by aweighting function ω₂({right arrow over (x)},t), and then processed bybackprojecting the differentiated projection data to obtain the functionQ({right arrow over (x)},{right arrow over (x)}′). In the second step,the function Q({right arrow over (x)},{right arrow over (x)}′) ismultiplied by a kernel function K({right arrow over (x)},{right arrowover (x)}′), and then integration is performed over the entire Euclideanspace spanned by the variable {right arrow over (x)}′, to obtain thereconstructed image value at the point {right arrow over (x)}.

It is important to elaborate on the meaning of the modified fan-beam andcone-beam projections g ^(A)[{circumflex over(β)}_({right arrow over (x)}′),{right arrow over (y)}(t)]. As shown inFIG. 8, the modified projections g ^(A)[{circumflex over(β)}_({right arrow over (x)}′),{right arrow over (y)}(t)] generallynon-zero in the whole Euclidean space. Therefore, one potential problemin the derived image reconstruction method is that we have tonumerically compute an integral over an infinite space. In addition,without an explicit choice of the weighting function, it is not yetclear whether the kernel function may posses an explicit and simpleform. Therefore, up to this point, it is not certain whether Eq.(2.24)-Eq. (2.26) provide a computationally efficient and practicalimage reconstruction method. Although we require the weighting functionto satisfy the factorization property (2.23), fortunately we still havesufficient degrees of freedom to choose an appropriate weightingfunction that will turn Eqs. (2.24)-(2.26) into a practical andefficient method. Namely, it is possible to obtain a one-dimensionalfiltering kernel. In this case, although the function Q({right arrowover (x)},{right arrow over (x)}′) is not band-limited along thefiltering lines, one can still use the data from a finite range toreconstruct the image. In the following paragraphs, we will show how tochoose the weighting function to make this kernel function simple andeasy to implement in practice. Before we discuss specific examples ofcone-beam and fan-beam image reconstruction via the FBPD, furtherclarification upon the properties of the weighting function isbeneficial. Since the weighting function satisfies the normalizationcondition (2.9), the components ω₁({right arrow over (x)},{circumflexover (k)}) and ω₂({right arrow over (x)},t) in Eq. (2.23) are notcompletely independent of one another. In fact, given a choice of thecomponent ω₂({right arrow over (x)},t), the component ω₁({right arrowover (x)},{circumflex over (k)}) is determined by the followingequation:

$\begin{matrix}{{\sum\limits_{t_{j} \in {T\mspace{11mu}{({\overset{->}{x},\hat{k}})}}}{{\omega_{1}\left( {\overset{->}{x},\hat{k}} \right)}\mspace{11mu}{\omega_{2}\left( {\overset{->}{x},t_{j}} \right)}\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot {{\overset{->}{y}}^{\prime}\left( t_{j} \right)}} \right\rbrack}}} = 1.} & (2.27)\end{matrix}$Since the function ω₁({right arrow over (x)},{circumflex over (k)}) doesnot depend on t, the summation over t_(j) for the last two terms onleft-hand-side of the above equation can be calculated by converting thesummation into an integral:

$\begin{matrix}{{{\sum\limits_{t_{j} \in {T\mspace{11mu}{({\overset{->}{x},\hat{k}})}}}{{\omega_{2}\left( {\overset{->}{x},t_{j}} \right)}\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot {{\overset{->}{y}}^{\prime}\left( t_{j} \right)}} \right\rbrack}}} = {{\int_{\Lambda}{{\mathbb{d}t}\;{\omega_{2}\left( {\overset{->}{x},t} \right)}\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot {{\overset{->}{y}}^{\prime}(t)}} \right\rbrack}{{\hat{k} \cdot {{\overset{->}{y}}^{\prime}(t)}}}\mspace{11mu}{\delta\mspace{11mu}\left\lbrack {\hat{k} \cdot \left( {\overset{->}{x} - {\overset{->}{y}(t)}} \right)} \right\rbrack}}} = {{\int_{\Lambda}{{\mathbb{d}t}{\int_{- \infty}^{+ \infty}\ {{\mathbb{d}l}\;{\omega_{2}\left( {\overset{->}{x},t} \right)}\mspace{11mu}{\hat{k} \cdot {{\overset{->}{y}}^{\prime}(t)}}\mspace{11mu}{\mathbb{e}}^{{\mathbb{i}2\pi}\; l\;{\hat{k} \cdot {\lbrack{\overset{->}{x} - {\overset{->}{y}{(t)}}}\rbrack}}}}}}} = {{\frac{1}{2{\pi\mathbb{i}}}\;{\int_{\Lambda}{{\mathbb{d}t}{\int_{- \infty}^{+ \infty}{\frac{dl}{l}\frac{\partial}{\partial t}\;{\mathbb{e}}^{{\mathbb{i}2\pi}\; l\;{\hat{k} \cdot {\lbrack{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}{(t)}}}\rbrack}}}{\omega_{2}\left( {\overset{->}{x},t} \right)}}}}}} = {\int_{\Lambda}{{\mathbb{d}t}\;{\omega_{2}\left( {\overset{->}{x},t} \right)}\frac{\partial}{\partial t}h\mspace{11mu}\left( {\overset{->}{x},\hat{k},t} \right)}}}}}},{where}} & (2.28) \\{{h\mspace{11mu}\left( {\overset{->}{x},\hat{k},t} \right)} = {{\frac{1}{2{\pi\mathbb{i}}}{\int_{- \infty}^{+ \infty}\ {{\mathbb{d}l}\frac{{\mathbb{e}}^{{\mathbb{i}2\pi}\; l\;{\hat{k} \cdot {\lbrack{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}{(t)}}}\rbrack}}}}{l}}}} = {\frac{1}{2}{{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot \left( {\overset{->}{x},{\overset{->}{y}\mspace{11mu}(t)}} \right)} \right\rbrack}.}}}} & (2.29)\end{matrix}$The function ω₁({right arrow over (x)},{circumflex over (k)}) is theinverse of the above integral value in Eq. (2.28). Namely, the factorω₁({right arrow over (x)},{circumflex over (k)}) is given by thefollowing equation:

$\begin{matrix}{{\omega_{1}^{- 1}\left( {\overset{->}{x},\hat{k}} \right)} = {\frac{1}{2}{\int_{\Lambda}\ {{\mathbb{d}t}\;{\omega_{2}\left( {\overset{->}{x},t} \right)}\frac{\partial}{\partial t}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot \left( {\overset{->}{x},{\overset{->}{y}\mspace{11mu}(t)}} \right)} \right\rbrack}}}}} & (2.30)\end{matrix}$Therefore, for a given image point {right arrow over (x)}, as long asone piece-wise constant function ω₂({right arrow over (x)},t) is chosenfor a segment of sufficient source trajectory, the factor ω₁({rightarrow over (x)},{circumflex over (k)}) is determined by the aboveequation. Consequently, the whole weighting function ω({right arrow over(x)},{circumflex over (k)}; t) is determined by equation (2.23). Afterthe factors ω₁({right arrow over (x)},{circumflex over (k)}) andω₁({right arrow over (x)},t) are specified, an image reconstructionalgorithm is given by equations (2.24)-(2.26).

We now develop an exact method for reconstructing an image from acquiredcone beam data. For the conventional helical source trajectory with 1-pidata acquisition, there is one and only one pi-line for each point inthe region-of-interest (ROI) defined by the cylindrical volume containedwithin the helix. However, in the case of an N-pi data acquisitionscheme, there are multiple different lines that pass through each imagepoint in the ROI and intersect the helical trajectory at two differentpoints. Recently, the pi-line concept has also been generalized to thecase of a helical trajectory with dynamical helical pitch and the caseof a helical trajectory with a variable radius. For these generalizedhelical trajectories, under certain conditions, there is one and onlyone pi-line for each point in the ROI inside a helix. In addition, thepi-line concept has been generalized to the saddle trajectory. Animportant feature of the saddle trajectory is that the pi-lines are notunique for a point within the ROI. However, a common feature of theaforementioned trajectories is the existence of at least one pi-line orone generalized pi-line. The pi-line provides a natural geometricguideline in choosing the weighting function. If we denote thecorresponding arc of a pi-line associated with an object point {rightarrow over (x)} as I({right arrow over (x)})=[t_(π) _(a) ({right arrowover (x)}), t_(π) _(b) ({right arrow over (x)})], we have a convenientchoice for the component ω₂({right arrow over (x)},t) in the weightingfunction:

$\begin{matrix}{{\omega_{2}\left( {\overset{->}{x},t} \right)} = \left\{ \begin{matrix}{1,{{{if}\mspace{14mu} t} \in {I\mspace{11mu}\left( \overset{->}{x} \right)}},} \\{0,{{otherwise}.}}\end{matrix} \right.} & (3.30)\end{matrix}$With this choice of ω₂({right arrow over (x)},t), the function ω₁({rightarrow over (x)},t) can be calculated as follows. Using Eq. (2.28), wehave

$\begin{matrix}{{\sum\limits_{t_{j} \in {T\mspace{11mu}{({\overset{->}{x},\hat{k}})}}}{{\omega_{2}\left( {\overset{->}{x},t_{j}} \right)}\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot {{\overset{->}{y}}^{\prime}\left( t_{j} \right)}} \right\rbrack}}} = {{\int_{I\mspace{11mu}{(\overset{->}{x})}}{{\mathbb{d}t}\frac{\partial}{\partial t}h\mspace{11mu}\left( {\overset{->}{x},\hat{k},t} \right)}} = {{{h\mspace{11mu}\left( {\overset{->}{x},\hat{k},t_{\pi_{a}}} \right)} - {h\mspace{11mu}\left( {\overset{->}{x},\hat{k},t_{\pi_{b}}} \right)}} = {\frac{1}{2}{\left\{ {{{sgn}\;\left\lbrack {\hat{k} \cdot \left( {\overset{->}{x} - {\overset{->}{y}\mspace{11mu}\left( t_{\pi_{a}} \right)}} \right)} \right\rbrack} - {{sgn}\left\lbrack {\hat{k} \cdot \left( {\overset{->}{x} - {\overset{->}{y}\left( t_{\pi_{b}} \right)}} \right)} \right\rbrack}} \right\}.}}}}} & (3.31)\end{matrix}$Since t_(π) _(a) and t_(π) _(b) are the endpoints of the pi-lineI({right arrow over (x)}), we havesgn[{circumflex over (k)}·({right arrow over (x)}−{right arrow over(y)}(t _(π) _(a) ))]−sgn[{circumflex over (k)}·({right arrow over(x)}−{right arrow over (y)}(t _(π) _(b) ))]=2sgn[{circumflex over(k)}·({right arrow over (y)}(t _(π) _(b) )−{right arrow over (y)}(t _(π)_(a) ))].  (3.32)Therefore, the function ω₁({right arrow over (x)},{circumflex over (k)})can be determined using Eq. (2.27):ω₁({right arrow over (x)},{circumflex over (k)})=sgn[{circumflex over(k)}·({right arrow over (y)}(t _(π) _(b) )−{right arrow over (y)}(t _(π)_(a) ))].  (3.33)

Using the explicit form of the component ω₁({right arrow over(x)},{circumflex over (k)}) in Eq. (3.33), the kernel function K({rightarrow over (x)},{right arrow over (x)}′) can be calculated from Eq.(2.25) as:

$\begin{matrix}{{{K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {\frac{1}{2{\pi\mathbb{i}}}\;{\int_{R_{k}^{3}}{{\mathbb{d}\overset{->}{k}}\mspace{11mu}{{sgn}\left\lbrack {\hat{k} \cdot {{\overset{->}{e}}_{\pi}\left( \overset{->}{x} \right)}} \right\rbrack}\;{\mathbb{e}}^{{\mathbb{i}2\pi}\;{\overset{->}{k} \cdot {({\overset{->}{x} - {\overset{->}{x}}^{\prime}})}}}}}}},} & (3.34)\end{matrix}$where {right arrow over (e)}_(π)({right arrow over (x)})={right arrowover (y)}(t_(π) _(b) )−{right arrow over (y)}(t_(π) _(a) ) has beenintroduced along the pi-line, we obtain:

$\begin{matrix}{{{{K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {\frac{1}{2\pi^{2}}\frac{1}{z_{\pi} = z_{\pi}^{\prime}}\delta\mspace{11mu}\left( {x_{\pi} - x_{\pi}^{\prime}} \right)\mspace{11mu}\delta\mspace{11mu}\left( {y_{\pi} - y_{\pi}^{\prime}} \right)}},}\mspace{14mu}} & (3.35)\end{matrix}$where {right arrow over (x)} and {right arrow over (x)}′ are now alongpi-lines, viz. {right arrow over (x)}=(x_(π),y_(π),z_(π)) and {rightarrow over (x)}′=(x′_(π),y′_(π),z′_(π)). Therefore, we obtain anexplicit one-dimensional filtering kernel. Here we show that the samekernel can be used in any source trajectory provided that the concept ofa pi-line is meaningful. Therefore, the filtering kernel is ashift-invariant Hilbert kernel along the pi-line. In fact, the same formof the function ω₂({right arrow over (x)},t) as given in Eq. (3.30) maybe chosen to reconstruct the points on an arbitrary straight line thatpasses one point and intersects the source trajectory at two differentpoints t_(a) and t_(b). As long as the form of ω₂({right arrow over(x)},t) is chosen as given in Eq. (3.30), the filtering kernel given inEq. (3.34) is applicable.

The backprojection function Q({right arrow over (x)},{right arrow over(x)}′) is given by:

$\begin{matrix}{{Q\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {{\int_{I(\overset{->}{x})}\ {{\mathbb{d}t}\frac{1}{{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}}\frac{\partial}{\partial q}{{\overset{\_}{g}}^{A}\left\lbrack {{\hat{\beta}}_{{\overset{->}{x}}^{\prime}},{\overset{->}{y}\mspace{11mu}(q)}} \right\rbrack}}}❘_{q = t}.}} & (3.36)\end{matrix}$The image point, {right arrow over (x)}, dependence of this functionshows up through the pi-line. The function g ^(A)[{circumflex over(β)}_({right arrow over (x)}′),{right arrow over (y)}(t)] is defined inEq. (2.21).

In summary, a mathematically exact image reconstruction method has beenderived for any source trajectory provided that the concept of a pi-lineor a generalized pi-line is meaningful and the source trajectorysatisfies the Tuy data sufficiency condition. In this method, an imageis reconstructed by filtering the backprojection image of differentiatedprojection data using a one-dimensional Hilbert kernel as given by Eq.(3.36) and Eq. (3.35).

Equation (3.30) gives one example of many possible selections ofweighting functions. A specific choice of the weighting functionω₂({right arrow over (x)},t) gives a specific image reconstructionformula. Namely, the filtering directions that are specified by theFourier transforms of the factors ω₁({right arrow over (x)},{circumflexover (k)}) as indicated in equation (2.25) may be different. In theabove example, the filtering direction is along a line that passesthrough the image point and intersects with the source trajectory at twodifferent points (a pi-line or a generalized pi-line). The weightingfactor ω₂({right arrow over (x)},t) may also be chosen in such a waythat the filtering direction is not along a pi-line. An example is givenas below:

$\begin{matrix}{{\omega_{2}\left( {\overset{->}{x},t} \right)} = \left\{ \begin{matrix}{1,{t \in \left\lbrack {t_{a},t_{0}} \right\rbrack}} \\{{- 1},{t \in \left\lbrack {t_{0},t_{b}} \right\rbrack}} \\{0,{otherwise}}\end{matrix} \right.} & (3.37)\end{matrix}$where t₀ is any point on the source trajectory. Using this weightingfactor ω₂({right arrow over (x)},t), the weighting factor ω₁({rightarrow over (x)},{circumflex over (k)}) is calculated as:

$\begin{matrix}{{{\omega_{1}\left( {\overset{->}{x},\hat{k}} \right)} = {{sgn}\mspace{11mu}\left\lbrack {{\hat{k} \cdot \hat{e}}\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}} \right)} \right\rbrack}},\mspace{14mu}{{\hat{e}\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}} \right)} = \frac{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}}{{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}}}}} & (3.38)\end{matrix}$

If we align the z-axis of the vector {right arrow over (k)} along theline labeled by the unit vector

$\;{{{\hat{e}\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}} \right)} = \frac{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}}{{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}\left( t_{0} \right)}}}},}$we obtain:

$\begin{matrix}{{K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {\frac{1}{2\pi^{2}}\frac{1}{z_{e} - z_{e}^{\prime}}\delta\mspace{11mu}\left( {x_{e} - x_{e}^{\prime}} \right)\mspace{11mu}\delta\mspace{11mu}\left( {y_{e} - y_{e}^{\prime}} \right)}} & (3.39)\end{matrix}$where {right arrow over (x)} and {right arrow over (x)}′ are now alongpi-lines, {right arrow over (x)}=(x_(e),y_(e),z_(e)) and {right arrowover (x )}′=(x′_(e),y′_(e),z′_(e)) backprojection image is given by:

$\begin{matrix}{{Q\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {{\left( {{\int_{t_{a}}^{t_{0}}{\mathbb{d}t}} - {\int_{t_{0}}^{t_{b}}{\mathbb{d}t}}} \right)\frac{1}{{\overset{->}{x} - {\overset{->}{y}\mspace{11mu}(t)}}}\frac{\partial}{\partial q}{{\overset{\_}{g}}^{A}\left\lbrack {{\hat{\beta}}_{\overset{->}{x}}^{\prime},(q)} \right\rbrack}}❘_{q = t}}} & (3.40)\end{matrix}$Note that the filtering direction is along a straight line that connectsimage point {right arrow over (x)} and pre-selected source point {rightarrow over (y)}(t₀). This filtering line only intersects the sourcetrajectory at one point {right arrow over (y)}(t₀). It is also importantto note that the selection of point {right arrow over (y)}(t₀) isarbitrary, therefore, we conclude that the filtering direction may beselected along any direction that passes through the image point {rightarrow over (x)}.

In summary, using the factorized weighting function in Eq. (2.23) and apiece-wise constant weighting factor ω₂({right arrow over (x)},t), animage reconstruction method may be developed using a one-dimensionalHilbert filtering of the backprojection image of differentiatedprojection data. This reconstruction method is applicable to any sourcetrajectory that fulfills the Tuy data sufficiency condition and api-line or a generalized pi-line is meaningful. The filtering directionmay be along any straight line that connects the image point and apreselected source point.

In this section, we consider image reconstruction from fan-beamprojections on an arc source trajectory with radius r. The sourcetrajectory is parameterized as:{right arrow over (y)}(t)=r(cos t, sin t),t∈[t _(i) , t _(f)],  (4.37)where t_(i) and t_(f) are view angles corresponding to the starting andending points on the scanning path. For this scanning path, there areinfinitely many lines that pass through an object point {right arrowover (x)} and intersect the source trajectory at two points. In order toderive an explicit FBPD reconstruction method, we discuss two possiblechoices of the component ω₂({right arrow over (x)},t) of the weightingfunction ω({right arrow over (x)},{circumflex over (k)}; t).

Case 1: Weighting function chosen for filtering along parallel lines

-   -   The intersection points of this horizontal line with the source        trajectory are labeled by {right arrow over (y)}(t_(a)({right        arrow over (x)})) and {right arrow over (y)}(t_(b)({right arrow        over (x)})). The first choice of the function ω₂({right arrow        over (x)},t) is:

$\begin{matrix}{{\omega_{2}\left( {\overset{->}{x},t} \right)} = \left\{ \begin{matrix}{1,{{{if}\mspace{14mu} t} \in \left\lbrack {t_{a},t_{b}} \right\rbrack},} \\{0,{{otherwise}.}}\end{matrix} \right.} & (4.38)\end{matrix}$

Using Eqs. (2.27) and (2.28), we obtain:ω₁({right arrow over (x)},{circumflex over (k)})=sgn[{circumflex over(k)}·({right arrow over (y)}(t _(b))−{right arrow over (y)}(t_(a)))].  (4.39)Substituting the above equation into Eq. (2.25), we obtain:

$\begin{matrix}{{K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {{- \frac{1}{2\pi^{2}}}\frac{1}{x - x^{\prime}}\delta\mspace{11mu}{\left( {y - y^{\prime}} \right).}}} & (4.40)\end{matrix}$Here we have aligned the vector {right arrow over (k)} along thehorizontal x-direction. Thus, the filtering is a Hilbert transformconducted along the horizontal direction. In the case of an equi-angulardetector, the derivative filtering backprojection function Q({rightarrow over (x)},{right arrow over (x)}′) is given by:

$\begin{matrix}\begin{matrix}{{{Q\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {{\int_{t_{a}(\overset{->}{x})}^{t_{b}(\overset{->}{x})}{{\mathbb{d}t}\frac{1}{{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}}\frac{\partial}{\partial q}{{\overset{\_}{g}}^{A}\left\lbrack {{\hat{\beta}}_{{\overset{->}{x}}^{\prime}},{\overset{->}{y}\mspace{11mu}(q)}} \right\rbrack}}}❘_{q = t}}},} \\{= {\int_{t_{a}(\overset{->}{x})}^{t_{b}(\overset{->}{x})}{{\mathbb{d}t}\frac{1}{{{\overset{->}{x}}^{\prime} - {\overset{->}{y}\mspace{11mu}(t)}}}\left( {\frac{\partial}{\partial t} - \frac{\partial}{\partial\gamma}} \right)}}} \\{\left\lbrack {{{gm}\mspace{11mu}\left( {\gamma,t} \right)}❘_{{\gamma = \gamma_{\overset{->}{x}}},}{{{- {gm}}\mspace{11mu}\left( {\gamma,t} \right)}❘_{{\gamma = \gamma_{\overset{->}{x}}},{+ \pi}}}} \right\rbrack,}\end{matrix} & (4.41)\end{matrix}$where g_(m)(γ,t) is the measured fan-beam projection data at view anglet and fan angle γ. The values of γ_({right arrow over (x)}), and{circumflex over (β)}_({right arrow over (x)}), are determined by:γ_({right arrow over (x)}′) =φ _({right arrow over (x)}′)−t+π,{circumflex over (β)} _({right arrow over (x)}′)=(cosφ_({right arrow over (x)}′), sin φ_({right arrow over (x)}′)),  (4.42)as demonstrated in FIG. 10.

The final reconstruction formula is given by:

$\begin{matrix}{{f\mspace{11mu}\left( \overset{->}{x} \right)} = {{- \frac{1}{2\pi^{2}}}{\int_{- \infty}^{+ \infty}\ {{\mathbb{d}x^{\prime}}{\int_{- \infty}^{+ \infty}\ {{\mathbb{d}y}\frac{1}{x - x^{\prime}}\delta\mspace{11mu}\left( {y - y^{\prime}} \right)\mspace{11mu} Q\mspace{11mu}{\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right).}}}}}}} & (4.43)\end{matrix}$For a given reconstruction point {right arrow over (x)}=(x,y), thebackprojection function Q({right arrow over (x)},{right arrow over(x)}′) should be calculated over an infinite horizontal line whichpasses through {right arrow over (x)}. Namely, {right arrow over(x)}′=(x′,y′) in Eq. (4.41) extends along an infinite one-dimensionalline. If the image function is compactly supported, a modified Hilberttransform may be used to perform the filtering process which onlyrequires data over a finite range.

The implementation steps are summarized below:

-   -   Calculate the backprojection image values of differentiated        fan-beam data on a Cartesian grid.    -   Filter the backprojection image with the Hilbert kernel along        horizontal lines as shown in FIG. 9.

The advantage of this method is that the image reconstruction can beconducted directly on a Cartesian grid. However, the disadvantage isthat, for different filtering lines, varying amounts of the projectiondata are used in the derivative backprojection step. Therefore, thischoice of weighting function potentially leads to nonuniform noisedistribution.

Case 2: Weighting function chosen for filtering along intersecting lines

-   -   In this case, the objective is to use the same amount of        projection data to reconstruct all the points within an ROI. For        a given image point {right arrow over (x)}, two lines are used        to connect the image point {right arrow over (x)} with the two        ending points at the view angles t_(i) and t_(f). As shown in        FIG. 11 these two straight lines intersect with the scanning        path at the view angles t_(a)({right arrow over (x)}) and        t_(b)({right arrow over (x)}) respectively. Consequently, the        scanning path is divided into three arcs:        T ₁({right arrow over (x)})={t|t _(i) <t<t _(a)({right arrow        over (x)})},        T ₂({right arrow over (x)})={t|t _(a)({right arrow over        (x)})<t<t _(b)({right arrow over (x)})},        T ₃({right arrow over (x)})={t|t _(b) <t<t _(f)}.  (4.44)        The function ω₂({right arrow over (x)},t) is chosen to be:

$\begin{matrix}{{\omega_{2}\left( {\overset{->}{x},t} \right)} = \left\{ \begin{matrix}{a,} & {{{{if}\mspace{14mu} t} \in {T_{1}\left( \overset{->}{x} \right)}},} \\{1,} & {{{{if}\mspace{14mu} t} \in {T_{2}\left( \overset{->}{x} \right)}},} \\{b,} & {{{{if}\mspace{14mu} t} \in {T_{3}\left( \overset{->}{x} \right)}},}\end{matrix} \right.} & (4.45)\end{matrix}$where constants a and be satisfy the condition:a+b=1, a≠b.  (4.46)

Using the above function ω₂({right arrow over (x)},t), the integral inthe Eq. (2.28) can be directly calculated to yield:

$\begin{matrix}{{\sum\limits_{t_{j} \in {T\mspace{11mu}{({\overset{->}{x},\hat{k}})}}}\;{{\omega_{2}\left( {\overset{->}{x},t_{j}} \right)}\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k} \cdot {{\overset{->}{y}}^{\prime}\left( t_{j} \right)}} \right\rbrack}}} = {{\int_{\Lambda}\ {{\mathbb{d}t}\;{\omega_{2}\left( {\overset{->}{x},t} \right)}\frac{\partial}{\partial t}h\mspace{11mu}\left( {\overset{->}{x},\hat{k},t} \right)}} = {{{- a}\mspace{11mu}{{sgn}\mspace{11mu}\left\lbrack {\hat{k}{\cdot \left( {{\overset{->}{y}\mspace{11mu}\left( t_{b} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{i} \right)}} \right)}} \right\rbrack}} + {b\mspace{11mu}{{sgn}\;\left\lbrack {{\hat{k} \cdot \left( {\overset{->}{y}\mspace{11mu}\left( {t_{a} - {\overset{->}{y}\mspace{11mu}\left( t_{f} \right)}} \right)} \right\rbrack},} \right.}}}}} & (4.47)\end{matrix}$where in the last step, the fact that the point {right arrow over (x)}lies at the intersection between the line passing through ({right arrowover (y)}(t_(i)), {right arrow over (y)}(t_(b))) and the line passingthrough ({right arrow over (y)}(t_(f)), {right arrow over (y)}(t_(a)))has been used.

The inverse of the above equation provides the factor ω₁({right arrowover (x)},{circumflex over (k)}) in the weighting function. The resultis:

$\begin{matrix}{{\omega_{1}\left( {\overset{->}{x},\hat{k}} \right)} = {{\frac{1}{a - b}{{sgn}\mspace{11mu}\left\lbrack {\hat{k}{\cdot \left( {{\overset{->}{y}\mspace{11mu}\left( t_{b} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{i} \right)}} \right)}} \right\rbrack}} + {\frac{b}{b - a}{{{sgn}\mspace{11mu}\left\lbrack {\hat{k}{\cdot \left( {{\overset{->}{y}\mspace{11mu}\left( t_{a} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{f} \right)}} \right)}} \right\rbrack}.}}}} & (4.48)\end{matrix}$

For simplicity, we introduce the following notation:

$\begin{matrix}{{{{\hat{e}}_{i}\left( \overset{->}{x} \right)} = \frac{{\overset{->}{y}\mspace{11mu}\left( t_{b} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{i} \right)}}{{{\overset{->}{y}\mspace{11mu}\left( t_{b} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{i} \right)}}}},{{{\hat{e}}_{f}\left( \overset{->}{x} \right)} = {\frac{{\overset{->}{y}\mspace{11mu}\left( t_{a} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{f} \right)}}{{{\overset{->}{y}\mspace{11mu}\left( t_{a} \right)} - {\overset{->}{y}\mspace{11mu}\left( t_{f} \right)}}}.}}} & (4.49)\end{matrix}$Substituting Eq. (4.48) into Eq. (2/25), we obtain,

$\begin{matrix}{{{K\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {{\frac{1}{2{\pi^{2}\left( {b - a} \right)}}\text{\{}\frac{a}{\left( {\overset{->}{x} - {\overset{->}{x}}^{\prime}} \right) \cdot {\hat{e}}_{i}}{\delta\mspace{11mu}\left\lbrack {\left( {\overset{->}{x} - {\overset{->}{x}}^{\prime}} \right) \cdot {\hat{e}}_{i}^{\bot}} \right\rbrack}} + {\frac{b}{\left( {\overset{->}{x} - {\overset{->}{x}}^{\prime}} \right) \cdot {\hat{e}}_{f}}{\delta\mspace{11mu}\left\lbrack {\left( {\overset{->}{x} - {\overset{->}{x}}^{\prime}} \right) \cdot {\hat{e}}_{f}^{\bot}} \right\rbrack}}}},} & (4.50)\end{matrix}$where {ê_(i) ^(⊥,ê) _(i)} and {ê_(f) ^(⊥),ê_(f)} are two sets oforthonormal bases,ê _(i) ^(⊥) ·ê _(i)=0, ê _(f) ^(⊥) ·ê _(f)=0.  (4.51)Therefore, the filtering kernel is still the Hilbert kernel. However,the filtering lines are along the two directions ê_(i) and ê_(f),respectively. In this case, the backprojection operation is given by:

$\begin{matrix}{{{Q\mspace{11mu}\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right)} = {\left\lbrack {a{\int_{T_{1{(\overset{->}{x})}}}{+ {\int_{T_{2{(\overset{->}{x})}}}{{+ b}\int_{T_{3{(\overset{->}{x})}}}}}}}} \right\rbrack\ {\mathbb{d}t}\ \frac{1}{{{\overset{->}{x}}^{\prime} = {\overset{->}{y}\mspace{11mu}(t)}}}{\left( {\frac{\partial}{\partial t} - \frac{\partial}{\partial\gamma}} \right)\left\lbrack {{g_{m}\left( {\gamma,t} \right)}❘_{\gamma = \gamma_{{\overset{->}{x}}^{\prime}}}{{- {g_{m}\left( {\gamma,t} \right)}}❘_{\gamma = \gamma_{{\overset{->}{x}}^{\prime} + x}}}} \right\rbrack}}},} & (4.52)\end{matrix}$where γ_({right arrow over (x)}), is determined by Eq. (4.42) and theweighting factors a and b satisfy Eq. (4.46).

The final reconstruction formula is given by:

$\begin{matrix}{{f\mspace{11mu}\left( \overset{->}{x} \right)} = {{- \frac{1}{2{\pi^{2}\left( {b - a} \right)}}}{\int_{- \infty}^{+ \infty}{{\mathbb{d}x^{\prime}}{\int_{- \infty}^{+ \infty}{{\mathbb{d}y}\left\{ {{a\frac{\delta\;\left\lbrack \left( {\overset{->}{x} - {{\overset{->}{x}}^{\prime} \cdot {\hat{e}}_{i}^{\bot}}} \right) \right\rbrack}{\left( {\overset{->}{x} - {\overset{->}{x}}^{\prime}} \right) \cdot {\hat{e}}_{i}}} + {b\frac{\delta\;\left\lbrack \left( {\overset{->}{x} - {{\overset{->}{x}}^{\prime} \cdot {\hat{e}}_{f}^{\bot}}} \right) \right\rbrack}{\left( {\overset{->}{x} - {\overset{->}{x}}^{\prime}} \right) \cdot {\hat{e}}_{f}}}} \right\}\mspace{11mu} Q\mspace{11mu}{\left( {\overset{->}{x},{\overset{->}{x}}^{\prime}} \right).}}}}}}} & (4.53)\end{matrix}$For a given reconstruction point {right arrow over (x)}=(x,y), thebackprojection function Q({right arrow over (x)},{right arrow over(x)}′) should be calculated over two infinite lines which pass through{right arrow over (x)}. Namely, {right arrow over (x)}′=(x′,y′) in Eq.(4.52) extends along two infinite one-dimensional lines, as shown inFIG. 11.

The implementation steps are summarized below:

-   -   Calculate the backprojection image values of differentiated        fan-beam data on a Cartesian grid.    -   Filter the backprojection image with the Hilbert kernel along        the lines labeled by ê_(i) and ê_(f) shown in FIG. 11.    -   Sum the weighted and filtered contributions from the two        filtering lines.

The advantage of this method is that the same amount of projection datais used for all the reconstructed points within an ROI. However, thedisadvantage is that, for a single point a filtering process has to beperformed along two filtering lines.

The weighting factors in Eqs. (3.37) and (3.38) may also be utilized inthe fan-beam case to develop fan-beam image reconstruction procedures.In this case, the filtering line is along the line connecting the imagepoint {right arrow over (x)} and a preselected source point {right arrowover (y)}(t₀).

Any other selections of a piece-wise constant weighting factor ω₂({rightarrow over (x)},t) will give a one-dimensional filtering kernel andthus, a new image reconstruction procedure by following the aboveexamples. Namely, from a given piece-wise constant weighting factorω₂({right arrow over (x)},t), we calculated the weighting factorω₁({right arrow over (x)},{circumflex over (k)}) using equation (2.20)and the backprojection image using equation (2.26). The filtering kernelis calculated according to Eq. (2.25) using the calculated weightingfactor ω₁({right arrow over (x)},{circumflex over (k)}).

The steps for numerical implementation of the above image reconstructionmethods are as follows:

Step 1: For a given direction {circumflex over(β)}_({right arrow over (x)}′), calculate the differentiation of theprojection data g ^(A)[{circumflex over(β)}_({right arrow over (x)}′),{right arrow over (y)}(t)] with respectto the source parameter t.Step 2: Calculate the backprojection image values Q({right arrow over(x)},{right arrow over (x)}′) for the points along the filteringdirection using Equation (2.26).Step 3: Numerically calculate the Hilbert transform of thebackprojection image along the filtering line to obtain image valuesƒ({right arrow over (x)}) using the following equation:

$\begin{matrix}{{f\mspace{11mu}\left( {x_{f},y_{f},z_{f}} \right)} = {{- \frac{1}{2\pi^{2}}}{\int_{- \infty}^{+ \infty}{{\mathbb{d}x_{f}^{\prime}}\frac{1}{x_{f} - x_{f}^{\prime}}Q\mspace{11mu}\left( {x_{f},y_{f},{z_{f};x_{f}^{\prime}},{y_{f}^{\prime} = y_{f}},{z_{f}^{\prime} = z_{f}}} \right)}}}} & (4.54)\end{matrix}$where the subscript “ƒ” means the Hilbert transform is one-dimensionaland is along the filtering direction.

Note that the backprojection image of differentiated projection data isnot compactly supported. Thus, a direct implementation of Hilbertfiltering in Equation (4.54) will introduce numerical errors.Fortunately, for a compactly supported image object, it is known thatsamples in a finite range may be utilized to accurately reconstruct theobject. Namely, the integral range in Equation (4.54) may be finite. Tobe more specific, the following modified formula can be utilized in anumerical implementation:

$\begin{matrix}{{f\mspace{11mu}\left( {x_{f},y_{f},z_{f}} \right)} = {\frac{1}{2\pi^{2}\sqrt{\left( {B - x_{f}} \right)\left( {x_{f} - A} \right)}}\left\lbrack \begin{matrix}{\int_{A}^{B}\ {{\mathbb{d}x_{f}^{\prime}}\frac{\sqrt{\left( {B - x_{f}^{\prime}} \right)\left( {x_{f}^{\prime} - A} \right)}}{x_{f} - x_{f}^{\prime}}}} \\{{Q\mspace{11mu}\left( {x_{f},y_{f},{z_{f};x_{f}^{\prime}},{y_{f}^{\prime} = y_{f}},{z_{f}^{\prime} = z_{f}}} \right)} + C}\end{matrix} \right\rbrack}} & (4.55)\end{matrix}$Where the constants A and B are integral limits which should be selectedoutside of the image support. The constant C is given by the projectiondata along the filtering direction from the following relation:C=−2πg _(m)({circumflex over (β)}_(f) ,t)|  (4.56)where {circumflex over (β)}_(f) is along the filtering direction.

When a flat-panel detector is utilized to detect the attenuated X-rays,the derivative may be calculated in terms of detector coordinates u, v,and source parameter t. Namely, the differentiated projection data isgiven by:

$\begin{matrix}{{g_{d}\left( {u,v,t} \right)} = {\left( {{\frac{D^{2} + u^{2}}{D}\frac{\partial}{\partial u}} + {\frac{uv}{D}\frac{\partial}{\partial v}} - \frac{\partial}{\partial t}} \right)\mspace{11mu}{{g_{m}\left( {u,v,t} \right)}.}}} & (4.57)\end{matrix}$If a colinear detector is utilized to detect the attenuated x-rays infan-beam case, then the index v is suppressed in equation (4.57).

When a curved CT-detector is utilized to detect the attenuated x-rays,the above derivative may also be calculated in terms of detectorcoordinates γ, α, and source parameter t. Where parameter γ is theparameter along the transverse direction (fan angle) and the parameter αis the cone angle.

$\begin{matrix}{{g_{d}\left( {\gamma,\alpha,t} \right)} = {\left( {\frac{\partial}{\partial\gamma} - \frac{\partial}{\partial t}} \right)\;{g_{m}\left( {\gamma,\alpha,t} \right)}}} & (4.58)\end{matrix}$If the detector row is one, this is just the conventional fan beam case.In this case, the cone angle index α is suppressed.

DESCRIPTION OF THE PREFERRED EMBODIMENT

With initial reference to FIGS. 6 and 7, a computed tomography (CT)imaging system 10 includes a gantry 12 representative of a “thirdgeneration” CT scanner. Gantry 12 has an x-ray source 13 that projects acone beam of x-rays 14 toward a detector array 16 on the opposite sideof the gantry. The detector array 16 is formed by a number of detectorelements 18 which together sense the projected x-rays that pass througha medical patient 15. Each detector element 18 produces an electricalsignal that represents the intensity of an impinging x-ray beam andhence the attenuation of the beam as it passes through the patient.During a scan to acquire x-ray projection data, the gantry 12 and thecomponents mounted thereon rotate about a center of rotation 19 locatedwithin the patient 15.

The rotation of the gantry and the operation of the x-ray source 13 aregoverned by a control mechanism 20 of the CT system. The controlmechanism 20 includes an x-ray controller 22 that provides power andtiming signals to the x-ray source 13 and a gantry motor controller 23that controls the rotational speed and position of the gantry 12. A dataacquisition system (DAS) 24 in the control mechanism 20 samples analogdata from detector elements 18 and converts the data to digital signalsfor subsequent processing. An image reconstructor 25, receives sampledand digitized x-ray data from the DAS 24 and performs high speed imagereconstruction according to the method of the present invention. Thereconstructed image is applied as an input to a computer 26 which storesthe image in a mass storage device 29.

The computer 26 also receives commands and scanning parameters from anoperator via console 30 that has a keyboard. An associated cathode raytube display 32 allows the operator to observe the reconstructed imageand other data from the computer 26. The operator supplied commands andparameters are used by the computer 26 to provide control signals andinformation to the DAS 24, the x-ray controller 22 and the gantry motorcontroller 23. In addition, computer 26 operates a table motorcontroller 34 which controls a motorized table 36 to position thepatient 15 in the gantry 12.

As shown best in FIG. 5, in a preferred embodiment of the presentinvention the detector array 16 is a flat array of detector elements 18,having N_(r) (e.g., 1000) elements 18 disposed along the in-plane (x,y)direction, and N_(z) (e.g., 16) elements 18 disposed along the z axis.The x-ray beam emanates from the x-ray source 13 and fans out as itpasses through the patient 15 and intercepts the detection array 16.Each acquired view is a N_(r) by N_(z) array of attenuation measurementsas seen when the gantry is oriented in one of its positions during thescan. The object of the present invention is to reconstruct a set of 2Dimage slices from the 3D array of acquired data produced by the x-raycone beam during the scan.

Because the cone beam diverges as it passes through the patient 15, thereconstruction of the parallel image slices is not possible with astraight forward fan beam filtering and backprojecting process. Thepresent invention enables an accurate reconstruction of the image slicesfrom this acquired cone beam data. In addition, the present inventionenables the accurate reconstruction of image slices even when the pathof the x-ray source 13 is less than an a complete circle or is a spiralpath or is a circle and a straight line path, or is two concentriccircles.

This reconstruction method is implemented in the image reconstructor 25.Referring particularly to FIG. 12, the cone beam projection data isreceived from the DAS 24 as a two-dimensional array of values which arepreprocessed in the standard manner at process block 40. Suchpreprocessing includes correcting for known errors and offsets andcalculating the minus log of the data to convert it to x-ray attenuationvalues.

As indicated at process block 42, the next step is to differentiate thepreprocessed attenuation profiles g_(m) (u,v,t) by calculating theweighted derivatives. With the flat detector shown in FIG. 5 this isdone according to the above Eq. 4.57, whereas Eq. 4.58 is employed whenthe attenuation data is acquired with a curved detector.

A loop is then entered in which each image pixel value f(x_(f)) iscalculated by a filtered backprojection process. As described above,this is done one filter line at a time, and the choice of the number anddirection of filter lines to be used in the image reconstruction willdepend on the particular scan performed and image prescribed. Morespecifically, the weighted backprojection Q(x_(f),x′_(f)) at each imagepixel along a filtering line is calculated as indicated at process block42. This is done according to Eq. (2.26) using projection views thatwere acquired over a segment of the source trajectory indicated by Eq.(??). The backprojected values are weighted by the product of weightingfactor ω₂({right arrow over (x)}′_(f),t) and the inverse of the distancebetween the backprojection image point and the x-ray source point(1/|{right arrow over (x)}′_(f)−{right arrow over (y)}(t)|) as indicatedin Eq. (2.26).

As indicated at process block 46, the backprojected values along thefilter line are filtered using a finite range Hilbert transform as setforth above in Eq. (4.55). This produces the final image values atpoints along the filter line and a test is then made at decision block48 to determine if additional image points need to be calculated. If so,the next filter line is selected as indicated at process block 50 andthe process repeats. As will be explained in more detail below with thespecific case studies, the nature of this test and the manner in whichthe next filter line is selected will depend on the particular strategyused to reconstruct an image from the acquired data set. Suffice to saythat the process repeats until the image values at locations throughoutimage space are calculated with sufficient density that an image ofprescribed resolution can be produced.

Depending on the particular reconstruction strategy chosen, thecalculated image values may or may not lie on a Cartesian grid. If not,as determined at decision block 52, a regridding process is used atprocess block 54 to calculate the image values at points on a Cartesiangrid. This is a well-known step used in medical imaging and it may be aninterpolation of the calculated image values or a convolution baseddistribution. The resulting Cartesian image may then be displayed andstored in Hounsfield units in the usual fashion as indicated at processblock 56.

Referring still to FIG. 12, in one preferred embodiment of the presentinvention the image reconstruction process is directed by a table 58which stores a set of filter line locations and associated weightingfactors ω₂({right arrow over (x)},t). These have been predetermined forthe prescribed scan as a list of filter lines and weighting functionsω₂({right arrow over (x)},t) which are properly supported by theacquired data set. The process described above repeats for each of thefilter lines listed in table 58 until all of the listed filter lineshave been used. This is determined at decision block 48. The directionas well as the location will be different for each filter line in thegeneral case.

There is a special case in which a series of filter lines can be definedwhich reconstruct the image as image points on a Cartesian grid. Thatis, the filter lines are all in the same direction, but are spaced apartfrom each other the same distance. This special case is for a fan beamacquisition as described above in case 1. In this embodiment theparallel filter lines are reconstructed one at a time until all of thelines needed to cover the prescribed image region of interest asdetermined at decision block 48.

Another preferred embodiment of the invention employs the presentinvention to reconstruct an image acquired with a fan beam configurationwith data truncation at every acquired projection view. This continuousview truncation (CVT) application is illustrated in FIG. 13 where anasymmetric detector 70 covers only one-half of the field of view of thefan beam produced by x-ray source 72. The detector 70 is asymmetric inthat the detector elements are located to one side of the iso-ray 76. Asthe x-ray source 72 and detector 70 revolve around the subject 74 in acircular path, only one-half the projection data is acquired at eachview angle.

Using the above-described image reconstruction method it proved possibleto accurately reconstruct an image of the entire subject 74 if thesource 72 and detector 70 performed a complete 27 revolution around thesubject 74. A series of radial filter lines were used that passed fromone boundary of the field of view and through the isocenter to theopposite boundary. In computing the backprojection image ofdifferentiated projection data, the derivatives were calculated by athree point formula. A Hann window was used to used to suppress the highfrequency fluctuations in the reconstructed images. The constant C inEq. 4.55 was determined by using Eq. 4.56. A bilinear interpolationscheme was used to find the projection values for a given filteringline. Namely, when the intersection of the filtering line and thecircular source trajectory is not right on a discrete focal spot, weneed to draw two parallel lines from two nearby focal spots and use theprojection values of these two parallel lines to interpolate theprojection value at the filtering line. However, these two parallellines may not hit on the detector cells. In order to obtain theprojection values at these two parallel lines, linear interpolation mayalso have to be introduced in the detector plane. Thus, we may haveinterpolation in view angles and in detector cells.

It has also been discovered that images of selected subregions of thesubject 74 can be accurately reconstructed from data sets acquired withCVT even when the scan covers less than a complete 2π revolution. Inthis case the filter lines are not radial, but a series of parallelfilter lines through the subregion at an angle that is supported byprojection views in the acquired data set. For example, if the subregionof interest is the heart which is a small part of the entire field ofview, the scan can be shortened to less than the short-scan requirementof 180° plus fan angle.

By imposing a factorization property on the weighting function ω({rightarrow over (x)},{circumflex over (k)};t) of a cone-beam inversionformula and its fan-beam counterpart, an image reconstruction method canbe defined as an integral of a product of two functions: K({right arrowover (x)},{right arrow over (x)}′) and Q({right arrow over (x)},{rightarrow over (x)}′). The function Q({right arrow over (x)},{right arrowover (x)}′) is the weighted backprojection of differentiated projectiondata acquired over a suitable range of view angles. The functionK({right arrow over (x)},{right arrow over (x)}′) is the Fouriertransform of one of the two weighting factors ω₁({right arrow over(x)},{circumflex over (k)}). By choosing the second weighting factorω₂({right arrow over (x)}, {circumflex over (k)}) to be piecewiseconstant, practical reconstruction methods result in which a filteringkernel is reduced to a Hilbert kernel or a linear combination of Hilbertkernels applied along different filtering lines. In order to calculatethe backprojection image Q({right arrow over (x)}_(f),{right arrow over(x)}′_(f)), only one x-ray projection from one view angle is required.Therefore, the detector size is not required to be large enough to coverthe whole image object or the whole transverse cross section and bothtransverse data truncation and longitudinal data truncations are allowedusing this reconstruction method. A detector that is not large enough tocover the transverse cross section of the image object or thelongitudinal cross section may be dynamically or adaptively translatedor moved during data acquisition to acquire attenuated projection dataprofiles for a complete backprojection image Q({right arrow over(x)}_(f),{right arrow over (x)}′_(f)).

A number of implementation variations are possible with the presentmethod. The differentiated projection data g_(d)(u,v;t) or g_(d)(γ,a;t)may be approximately calculated using the numerical difference betweentwo nearby data points. The differentiated projection data g_(d)(u,v;t)or g_(d)(γ,a;t) may also be approximately calculated using fast Fouriertransforms (FFT) in the Fourier domain. The Hilbert transform may becalculated in the real space domain using a shifted and addition method,or the Finite Hilbert transform may be calculated in the Fourier domainusing FFT.

This image reconstruction method may be utilized in other tomographicimaging modalities such as MRI, PET, SPECT, thermalacoustic CT (TCT),and photoacoustic CT (PCT). The detectors may be flat panel CTdetectors, curved CT detectors, gamma cameras (SPECT) or a transducer(TCT and PCT), or read-out coils (MRI). The constant C in the finiterange Hilbert transform may also be determined using the physicalinformation that the image values should be zero outside the imagesupport. Namely, when an image point is selected outside the imagesupport, the constant C is given by the negative value of the integralin the square bracket of Eq. (4.55). Also, the integral limit in thefinite Hilbert transform Equation (4.55) may be selected to be symmetric(A=−B).

It should be apparent that the selection of filtering line location andorientation is flexible. For a given targeted reconstruction region ofinterest or volume, one may select different filtering lines to fillthat entire image volume or region. The filtering line direction may beselected along a pi-line, a generalized pi-line, or an arbitrary linethat connects the image points and a preselected source position (openchord). For any image point inside a volume or region of interest, thereshould exist at least one straight line that passes through the imagepoint and intersects the source trajectory at two points in its scantrajectory. Otherwise, the method is not limited to particular scantrajectories or detector geometries.

1. A method for producing an image with a tomographic imaging system,the steps comprising: a) performing a scan in which a plurality of beamsare produced by a source and corresponding projection views of a subjectare acquired at different view angles; b) differentiating the acquiredprojection views; c) producing image values by backprojecting adifferentiated projection view along a filter line that passes throughimage space and weighting the backprojected value by the product of aweighting factor and the inverse of the distance between the projectionview source and the image value location; d) filtering the image valuesalong the filter line; and e) repeating steps c) and d) to produceadditional filtered image values along different filter lines throughouta region of interest in image space.
 2. The method as recited in claim 1which step d) employs a Hilbert transformation over a finite range. 3.The method as recited in claim 1 in which the filter lines aresubstantially parallel to each other.
 4. The method as recited in claim3 in which the filtered image values along the filter lines aresubstantially uniformly distributed throughout the region of interest.5. The method as recited in claim 1 which includes: f) regridding thefiltered image values such that they are substantially uniformlydistributed throughout the region of interest.
 6. The method as recitedin claim 1 in which the filter lines are radially directed and passthrough the center of image space at different angles.
 7. The method asrecited in claim 6 which includes: f) regridding the filtered imagevalues such that they lie on the Cartesian grid in image space.
 8. Themethod as recited in claim 1 in which step a) is performed by detectinga divergent beam.
 9. The method as recited in claim 8 in which thedivergent beam is a fan beam.
 10. The method as recited in claim 8 inwhich the divergent beam is a cone beam.
 11. The method as recited inclam 8 in which the divergent beam is detected symmetrically about aniso-ray.
 12. The method as recited in claim 8 in which the divergentbeam is detected asymmetrically about an iso-ray.
 13. The method asrecited in claim 8 in which the divergent beam is produced by an x-raysource moving in a path around the subject.
 14. A method for producingan image with a tomographic imaging system, the steps comprising: a)performing a scan in which a plurality of divergent beam projectionviews of a subject are acquired at different view angles; b)differentiating the acquired projection views; and c) producing an imagefrom the acquired projection views which is an integral of a product oftwo functions, K and Q, where the function Q is a weightedbackprojection operation on the differentiated projection view data,where the function K is a Fourier transform of a first weighting factorω₁, that is one of two weighting factors ω₁ and ω₂ that result fromfactorizing a weighting function ω.
 15. The method as recited in claim14 in which the weighting function ω₂ is piecewise constant.